from scipy import *
#It is easy to not use the quad
from scipy.integrate import quad
from numpy import *

k = 50*pi

# contour parametrization and its derivatives [alpha,betta]->R
alpha = 0.0
betta = 1.0
def ro1(t):    return sin(t)
def z1(t):    return cos(t)
def dro1(t):    return cos(t)
def dz1(t):    return -sin(t)
def ddro1(t):    return -sin(t)
def ddz1(t):    return -cos(t)
def dddro1(t):    return -cos(t)
def dddz1(t):    return sin(t)

#discretization order
length = quad(lambda t: sqrt(dro1(t)**2+dz1(t)**2),alpha,betta)
n = int(3*length[0]*k/pi) + 5.
print "n =",n

# the  roots of Chebyshev polynomials(Chebyshev points)
ind1 = arange(int(n-1))
ind2 = arange(int(n))
t0 = cos(pi*(ind1+1)/n)
t = cos(pi*(2*ind2+1)/(2*n))
ind = arange(2*n-1)
t0t = hstack((t0,t))

# values of contour parametrization 
# and its derivatives [-1,1]->R in the Chebyshev points
mult = (betta-alpha)/2
tau = alpha*(1-t0t)/2 + betta*(1+t0t)/2
ro = ro1(tau)
z = z1(tau)
dz = dz1(tau)*mult
dro = dro1(tau)*mult
ddz = ddz1(tau)*mult**2
ddro = ddro1(tau)*(mult**2)
dddro = dddro1(tau)*(mult**3)
dddz = dddz1(tau)*(mult**3)

#Lyame koefficient
h_tau = sqrt(dro**2+dz**2)

#distance between to point on the surface
def L(m,p,psi): return sqrt(ro[p]**2+ro[m]**2-2*ro[p]*ro[m]*cos(psi)+(z[p]-z[m])**2)

#FUCKING LIMITS

def f_Smooth_Re(M,p,psi):
    if psi==0 : return 0   
    else:
        LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))
        Part11 = -(dro[p]**2)*cos(psi)-dz[p]**2
        Part12 = -(cos(k*LL)/(LL**3) + k*sin(k*LL)/(LL**2))*cos(M*psi)
        Part13 = (1/LL**3)*(1 - (M**2 * sin(psi)**2)/2)
        Part14 = (k**2/2 + M**2/(2*ro[p]**2)*(cos(psi) - 1))/LL
        Part21 = (ro[p]*dro[p]*(1 - cos(psi)))**2
        Part22 = (cos(k*LL)/(LL**3)*(-k**2 + 3/(LL**2)) + 3*k*sin(k*LL)/(LL**4))*cos(M*psi)
        Part23 = -3*(1 - (M**2 * sin(psi)**2)/2)/(LL**5)
        Part24 = -(1/LL**3)*(k**2/2 + M**2/(2*ro[p]**2)*(cos(psi)-1))
        Res = Part11*(Part12 + Part13 + Part14) + Part21*(Part22+Part23+Part24)
        return Res
   
def f_Smooth_Im(M,p,psi):
    if psi==0.0 : return -k**3*(dro[p]**2+dz[p]**2)/3
    else: 
       LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))
       Part11 = -(dro[p]**2)*cos(psi)-dz[p]**2
       Part12 = (sin(k*LL)/(LL**3) - k*cos(k*LL)/(LL**2))*cos(M*psi)
       Part21 = (ro[p]*dro[p]*(1 - cos(psi)))**2
       Part22 = (-sin(k*LL)/(LL**3)*(-k**2 + 3/(LL**2)) + 3*k*cos(k*LL)/(LL**4))*cos(M*psi)
       Res = Part11*Part12 + Part21*Part22
       return Res

def Smooth(M,p):   
    Int_Re,err = quad(lambda psi: f_Smooth_Re(M,p,psi),0,pi)
    Int_Im,err = quad(lambda psi: f_Smooth_Im(M,p,psi),0,pi)
    return 2*Int_Re + 2*Int_Im*1j

def Lim_1_L(p):
    Part1 = (-dro[p]**2-dz[p]**2)/(2*ro[p]**2)
    mult1 = ddro[p]**2/2 + dro[p]*dddro[p]/3 + ddz[p]**2/2 + dz[p]*dddz[p]/3
    mult2 = dro[p]**2 + dz[p]**2
    mult3 = (dro[p]*ddro[p] + dz[p]*ddz[p])**2
    Part2 = -(mult1*mult2-mult3)/mult2**2
    Part3 = -dro[p]/ro[p]**2*(dro[p]/ro[p]-(dro[p]*ddro[p]+dz[p]*ddz[p])/mult2)
    Part4 = (dz[p]**2 + 2*dro[p]**2)/(2*ro[p]**3)*log(8*ro[p]/sqrt(mult2))    
    Part5 = mult2/(2*ro[p])**3*(5-2*log(8*ro[p]/sqrt(mult2)))   
    Part6 = -(3*dz[p]**2 + 2*dro[p]**2 + 2*ro[p]*ddro[p])/(4*ro[p]**3)   
    return (Part1 + Part2)/ro[p] + Part3 + Part4 + Part5+ Part6

def Lim_L(p):
    mult = dro[p]**2 + dz[p]**2    
    return -2*dz[p]**2/ro[p] + 2*mult/ro[p]*(2 - log(8*ro[p]/sqrt(mult)))

def Lim_ddS(M,p): 
    return Smooth(M,p) + Lim_1_L(p) + (M**2/(2*ro[p]**2) - k**2/2)*Lim_L(p)

def f_A_Re(M,p,psi):
    LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))    
    return (cos(k*LL)*cos(M*psi)-1)/LL

def f_A_Im(M,p,psi):
    LL = sqrt(2)*ro[p]*sqrt(1-cos(psi))    
    return -sin(k*LL)*cos(M*psi)/LL

def Smooth2(M,p):   
    Int_Re,err = quad(lambda psi: f_A_Re(M,p,psi),0,pi)
    Int_Im,err = quad(lambda psi: f_A_Im(M,p,psi),0,pi)
    return 2*Int_Re + 2*Int_Im*1j

def Lim_S(M,p): 
    return Smooth2(M,p) + log(8*ro[p]/sqrt(dro[p]**2 + dz[p]**2))*2/ro[p]

def Lim_K1tau(M,p):
    Part = dro[p]**2*(Lim_S(M+1,p)+Lim_S(M-1,p))/2 + dz[p]**2*Lim_S(M,p)    
    return ro[p]**3*(Lim_ddS(M,p) - k**2*Part)

#print Lim_K1tau(1,3)
#quadrature formulas for the hypersingular, singular integrals and integrals with the logarithm kernels

def delta(x, y):  return int(x == y)

v_delta = vectorize(delta)

def A(l,j): 
    if l == j: return -n/2
    else: return (1-(-1)**(j+l))*(1-t0[j]**2)/(t0[l]-t0[j])**2/n

def Gamma1(l,j):
    if l == j: return 0
    else: return (1-(-1)**(j+l))*(1-t0[j]**2)/(t0[l]-t0[j])/n

def T(p,t): 
    return cos(p*arccos(t))

def LI(tt,j):
    part = (T(ind1+1,t0[j])*T(ind1+1,tt)/(ind1+1)).sum()
    return -1/n*(1-t0[j]**2)*(log(2) + 2*part + (-1)**j*T(n,tt)/n)

def LII(tt,j):
    part = (T(ind1+1,t[j])*T(ind1+1,tt)/(ind1+1)).sum()
    return -1/n*(log(2) + 2*part)

def Gamma2(l,j): 
    return (1/n)*1/(t0[l]-t[j])

def Gamma3(l,j): 
    return (1/n)*(1-t0[j]**2)/(t0[l]-t[j])

#quadrature formulas for smooth kernels
#simplify
def m0(M,p,j):
    Lmin = sqrt((ro[j]-ro[p])**2 + (z[j]-z[p])**2)
    Lmax = sqrt((ro[j]+ro[p])**2 + (z[j]-z[p])**2)
    a = (Lmax - Lmin)/(3*Lmin + Lmax)   
    arg = 10**(-4)*(1-a**2)/pi            
    nn = abs(log10(arg)/log10(a)) + 2
    return round((Lmax-Lmin)*k/pi + nn + M)

#print m0(1,yy,xx)
def a1_tau(p): 
    return -2*ro[p]**2

def b1_tau(p): 
    return -dro[p]*ro[p]

def c1_tau(M,p): 
    return (k**2*(ro[p])**2 + M**2)*(dro[p]**2 + dz[p]**2) - (dz[p]**2 + 3*dro[p]**2)/4

def dg(p,j,psi): 
    return ro[j]*dro[j] - ro[p]*dro[j]*cos(psi) + (z[j] - z[p])*dz[j]

def d2g(p,j,psi):
    return -dro[p]*dro[j]*cos(psi) - dz[p]*dz[j]

def f_A_1tau(M,p,j,psi):
    LL = L(p,j,psi)
    Part1 = d2g(p,j,psi)*(1/LL**2)*(-1/LL - 1j*k)
    Part2 = dg(p,j,psi)*dg(j,p,psi)*1/LL**3*(-k**2 + 3*1j*k/LL + 3/LL**2)
    Part3 = -k**2*(dro[p]*dro[j]*cos(psi)/LL + dz[p]*dz[j]/LL)
    return exp(-1j*k*LL)*(Part1 + Part2 + Part3)*cos(M*psi)


#PROBLEM - IT IS WRONG ARG ARANGE(m), if m IS VECTOR


def A1_tau(M,p,j):
    m = int(m0(M,p,j))
    phi = 2*pi*(arange(m)+1)/(2*m+1)      
    part = sum(f_A_1tau(M,p,j,phi))   
    return 2*pi/(2*m+1)*(f_A_1tau(M,p,j,0) + 2*part)

def K1_tau(M,p,j): 
    if p!=j:    
        Part1 = ro[p]**3*A1_tau(M,p,j) - a1_tau(p)/(t0[p]-t0[j])**2 
        Part2 = - b1_tau(p)/(t0[p]-t0[j]) - c1_tau(M,p)*log(abs(t0[p]-t0[j]))       
        return (Part1 + Part2)*(1-t0[j]**2)/n
    else:
        return Lim_K1tau(M,p)*(1-t0[p]**2)/n
        
a = zeros((n-1,n-1),complex)


from time import time
t1 = time()
for p in xrange(a.shape[0]):
    for j in xrange(a.shape[1]):
        if j!=p :
            a[p,j]=K1_tau(1,p,j)
t2 = time()
print t2-t1

#print round_(a)

d = input()


